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We discuss that hadron-induced atmospheric air showers from ultra-high energy 
cosmic rays are sensitive to QCD interactions at very small momentum fractions 
x where nonlinear effects should become important. The leading partons from 
the projectile acquire large random transverse momenta as they pass through the 
strong field of the target nucleus, which breaks up their coherence. This leads to 
a steeper x^-distribution of leading hadrons as compared to low energy collisions, 
which in turn reduces the position of the shower maximum X max . We argue that 
high-energy hadronic interaction models should account for this effect, caused by 
the approach to the black-body limit, which may shift fits of the composition 
of the cosmic ray spectrum near the GZK cutoff towards lighter elements. We 
further show that present data on X max (E) exclude that the rapid ~ l/rr. 0,3 growth 
of the saturation boundary (which is compatible with RH1C and HERA data) 
persists up to GZK cutoff energies. Measurements of pA collisions at LHC could 
further test the small-x regime and advance our understanding of high density 
QCD significantly. 



1. Introduction 

Today, quite little is known about the origin, the spectrum and the composi- 
tion of the highest energy cosmic rays. For example, AGASA 1 found about 
10 events with E > 10 11 GeV, well above the Greisen-Zatsepin-Kuzmin 
(GZK) cutoff, Eqzk — 6 • 10 10 GeV, which arises because of interaction 
of protons with the cosmic microwave background. On the other hand, 
the results of the HIRES 2 collaboration agree with the existence of the 
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GZK cutoff (assuming isotropic sources). Forthcoming Auger 3 data near 
the GZK cutoff will provide higher statistics and hopefully help to resolve 
this puzzle. 

The precise knowledge of the primary cosmic ray properties, that is the 
particle type, energy and arrival direction, is crucial for the interpretation 
of their possible source and acceleration mechanism. Standard candidates 
for the highest energy cosmic rays are protons or heavier nuclei, being 
accelerated in extreme astrophysical phenomena, or photons, arising for 
example from the decay of ultra-heavy X— particles. 

Experiments detect cosmic rays indirectly via air showers induced when 
they enter the atmosphere. One tries to deduce the properties of the pri- 
mary particle from those of the induced shower. Therefore, a good un- 
derstanding of the physics of high-energy interactions in the atmosphere 
is mandatory. However, the maximum energies exceed those of terrestrial 
accelerators by far, and so our knowledge of hadronic interactions needs 
to be extrapolated to unknown regimes. Also, as will be discussed in more 
detail below, air showers are mostly sensitive to forward particle production 
which is less well measured in accelerator experiments. 

Several features of strong interactions are expected to change dramat- 
ically at very high energies. First, the parameters of the soft interactions 
change - the total cross section changes by a factor of <~ 3, while average im- 
pact parameters increase by ~ 50%. The changes in the (semi-)hard interac- 
tions are even more dramatic. Indeed, a leading parton from the projectile 
propagates through the very strong gluon fields in the target. For example, 
for a "low-energy" p + A collision with i?Lab = 400 GeV, a parton from the 
projectile carrying a momentum fraction x p ~ 0.1 receiving a transverse 
kick of p t ~ 2 GeV interacts with a gluon with xa = ^p\jx v s ~ 0.1. At 
GZK energy, -ELab ~ 10 11 GeV, this corresponds to xa ~ 10 -10 while direct 
measurements at HERA covered only the range x > 10~ 3 and even indirect 
ones are sensitive only down to x > 10~ 4 . This is six orders of magnitude 
above the x-range to which cosmic rays near the cutoff are sensitive. 

Studies at HERA indicate that the gluon density of the nucleon, 
xgj\[(x,Q 2 ), grows very strongly with decreasing momentum fraction x, 
roughly as x~ x( - q2 \ with A(Q 2 ) > 0.2 for Q 2 > 2 GeV 2 . The data can be 
fitted by the NLO QCD evolution equations. The analysis of partial waves 
for the interaction of a small dipole with the nucleon at HERA energies 
indicates that for qq dipoles the partial waves remain substantially below 
the unitarity limit, while for gg dipoles the unitarity limit is practically 
reached for virtualities Q 2 ~ 4 GeV 2 at top HERA energies. This indicates 
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that for higher energies further growth of the gluon density in the proton 
should be tamed for a range of virtualities that increases with energy. 

The rate at which the growth of the gluon density is "tamed" , depends 
strongly on its behavior at small x. One of the most popular approaches 
for a long time was the BFKL approximation where one sums the series of 
leading, next to leading log(l/x). The series was found to converge poorly, 
to a large extent due to the specifics of the treatment of the energy conser- 
vation effects. More recent calculations 4 ' 5 (some of which were discussed 
at the Erice meeting) try to treat simultaneously logarithms of both 1/x 
and Q 2 and to treat more accurately the phase space available for gluon 
emission at a given energy. They appear to indicate that the NLO DGLAP 
approximation should effectively work for x > 10 -3 for the scattering off 
a gluon or, correspondingly, for x > 1CP 4 for the scattering off a nucleon, 
which is consistent with the HERA findings. 

It appears natural to expect that the taming effects would still allow 
the interactions to reach the maximal possible strength allowed by unitarity 
over a wide range of impact parameters which should increase with energy. 
Indeed, this is implemented in all models currently on the market with the 
only difference being the rate of the approach to the unitarity limit. 

If the energies are large enough, the constituents of the projectile 
hadron/photon propagating through the nucleus resolve strong small- a; 
gluon fields in the target. During the propagation through such media 
they should experience strong distortions - at the very least they should 
obtain significant transverse momenta inversely proportional to the size of 
test dipoles for which the interaction becomes black. Also, some of the 
processes relevant in this case, like hard scattering of the projectile partons 
off small- a; partons, lead to fractional energy losses. 

However the most important effect for the purposes of near-GZK in- 
elastic collisions is loss of coherence of the leading partons of the projectile 
as they acquire random transverse momenta. This leads to independent 
fragmentation of the leading partons from the projectile over a large range 
of rapidities, and hence to a much softer energy spectrum of the leading 
particles 6,7 . 

In this paper we review our first efforts to model this effect and to 
show its relevance for the understanding of air showers induced by cosmic 
rays near the GZK cutoff. Our primary goal is to analyze the implication 
of various models for the small- a: behavior of the gluon densities beyond 
the HERA range. We demonstrate that already the current data on the 
longitudinal and lateral structure of giant air showers allow us to rule out 



February 2, 2008 2:1 Proceedings Trim Size: 9in x 6in 



ericcQCD 



4 

certain models where gluon densities increase very rapidly with energy 8 . 
At the same time, models which appear to be consistent with the recent 
theoretical studies 4,5 lead to relatively small effects which are consistent 
with the air shower data and suggest that the spectrum near the cutoff is 
dominated by protons. 

Finally, we also point out that the gluon densities encountered in central 
proton-nucleus collisions at LHC are similar to those for central proton-air 
collisions near the cutoff energy Hence, we also present some predictions for 
p-" heavy nucleus" central collisions at RHIC and LHC energies which could 
test our suggestion for the mechanism of energy degradation by projectile 
breakup. Such measurements could help us understand the interactions of 
very high energy cosmic rays with the atmosphere and, consequently, of 
their composition and origin. 

2. Scattering at high energies 

The wave function of a hadron (or nucleus) boosted to large rapidity ex- 
hibits a large number of gluons at small x. The density of gluons per unit of 
transverse area and of rapidity at saturation is denoted by Q 2 , the so-called 
saturation momentum. This provides an intrinsic momentum scale which 
grows with atomic number (for nuclei) and with rapidity, due to continued 
gluon radiation as phase space grows. For sufficiently high energies and/or 
large nuclei, Q s can grow much larger than Aqcd and so weak coupling 
methods are applicable. Nevertheless, the well known leading-twist pQCD 
can not be used when the gluon density is large; rather, scattering ampli- 
tudes have to be resummed to all orders in the density. When probed at 
a scale below Q s , scattering cross sections approach the geometrical size 
of the hadron (the "black body" limit). A perturbative QCD based mech- 
anism for unitarization of cross sections is provided by gluon saturation 
effects 9 > 10 > 1:L . On the other hand, for Q 2 ^> Q 2 the process occurs in the 
dilute DGLAP 12 regime where cross sections are approximately determined 
by the known leading- twist pQCD expressions. 

In this section we discuss particle production in the collision of a hadron, 
which for the present purposes is either a nucleon or a meson, with a target 
nucleus in the atmosphere. The development of the air shower is sensitive 
mainly to the distribution of the most energetic particles (see section 4) 
while the low-energy particles produced near the nuclear fragmentation 
region are less important for the observables studied here. Due to QCD 
evolution (section 2.3) this so-called forward region probes the high gluon 
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density (small-x) regime of the target, while the density of the projectile 
is rather low. Hence, in the relevant rapidity (or Feynman-Xi?) region we 
are dealing with a "dilute" projectile hadron impinging on a "dense" target 
nucleus: Q h s < Qf. 

This assumption breaks down at large impact parameters, where even 
the saturation momentum of the nucleus, Qf(y,b), is not large. For such 
events, as well as for collisions at moderately high energies, no intrinsic 
semi-hard scale exists in the problem and so a treatment within weak cou- 
pling QCD is not applicable. In accelerator experiments this regime could 
be avoided by appropriately tuning the control parameters, such as collision 
energy, atomic number of projectile and target, impact parameter (trigger), 
rapidity y, transverse momentum p t and so on. This, of course, is not fea- 
sible in the case of cosmic ray air showers; here, we model such collisions 
using the SIBYLL leading- twist event generator. This is discussed in more 
detail in section 3. 

2.1. Leading quarks 

With this in mind, we now focus on particle production in collisions at suffi- 
ciently high energy and sufficiently small impact parameter where the satu- 
ration momentum of the nucleus is large enough to warrant a weak-coupling 
approach. The dominant process for fast particle production (xp > 0.1) is 
scattering of quarks from the incident dilute projectile on the dense target. 
For high quark energy we assume that the eikonal approximation applies 
such that p + is conserved. The transverse momentum distribution of scat- 
tered quarks is then given by the correlation function of two Wilson lines, 



running up and down the light cone at transverse separation r t (in the 
amplitude and its complex conjugate), 



Here, the convention is that the incident hadron has positive rapidity, i.e. 
the large component of its light-cone momentum is P + , and that of the 
incoming quark is p + = xP + (q + for the outgoing quark). The two-point 
function has to be evaluated in the background field of the target nucleus. 
When this field is weak, the Wilson lines can be expanded to leading order 




(1) 



-qA _ 
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and the problem reduces to evaluating the two-point function of the gauge 
field A*. 

In the strong-field regime gA + ~ 1, however, one needs to evaluate 
the correlation function to all orders (corresponding to summing over any 
number of scatterings of the incident quark). A relatively simple closed 
expression can be obtained in the McLerran-Venugopalan model of the 
small- x gluon distribution of the dense target 10 . In that model, the small- a; 
gluons are described as a stochastic classical non-abelian Yang-Mills field 
which is averaged over with a Gaussian distribution, n-point functions then 
factorize into powers of the two-point function. The qA cross section is then 
given by 13 

■»+ 



Q 



dq+d 2 q t d 2 b 
C{qt) 



Q 
P+ 



P 



C{qt) 



(2tt) 2 

— 2 exp 



exp 



-Ql 
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n 



1 



(3) 



(4) 



This expression is valid to leading order in a s (tree level) , but to all orders 
in Q s since it resums any number of scatterings of the quark in the strong 
field of the nucleus. The saturation momentum Q s , as introduced in cq. (3), 
is related to x, the total color charge density squared (per unit area) from 
the nucleus integrated up to the rapidity y of the probe (i.e. the projectile 
quark), by Q 2 S = A-K 2 a 2 s \ (^c ~ 1)/N C . In the low-density limit, \ is 
proportional to the ordinary leading-twist gluon distribution function of 
the nucleus 14 : 

A 



x( X ) = ^ [ d*> Ql) + j^j 9( X ', Ql)) , 



(5) 



where q{x) and g(x) denote the quark and gluon distributions of a nucleon, 
respectively; note that shadowing in the linear regime 15 would tend to 
reduce x somewhat but is neglected here since we are dealing with small 
nuclei (mass number 14 — > 16) and because the induced air showers are 
sensitive mainly to the small- a; regime in the nucleus. 

The integrals over p t in eq. (3) are cut off in the infrared by some cutoff 
A, which we assume is of order Aqcd- At large transverse momentum, 
again the first exponential in (3) can be expanded order by order 16 ' 17 to 
generate the usual power series in 1/qf: 



1 

2V2 q f 



7T q 2 A 



C l 2 
Qt 



(6) 
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This expression is valid to leading logarithmic accuracy. The first term 
corresponds to the perturbative one-gluon t-channel exchange contribution 
to qg — > qg scattering and exhibits the well-known power-law divergence of 
leading-twist perturbation theory for small momentum transfer. 

On the other hand, for Q s > q t one obtains in the leading logarithmic 
approximation 17 



This approximation reproduces the behavior of the full expression (3) about 
qt ~ Qs, and hence the transverse momentum integrated cross section rea- 
sonably well. It is useful when the cutoff A <C Q s , that is, when color 
neutrality is enforced on distance scales of order 1/A > 1/Q S . If, how- 
ever, color neutrality in the target nucleus occurs over distances of order 
1/Qs 18 then A ~ Q s and one has to go beyond the leading-logarithmic 
approximation. 

It is essential to realize that the high-energy part of the air shower is 
essentially one-dimensional, i.e. the transverse momenta of the produced 
hadrons play no role a (see section 4). This, in turn, implies that when Q s 
is large that the high-transverse momentum leading-twist regime can be 
neglected. The ^-distribution of forward valence quarks can thus be taken 
to be given by the simple expression (7) rather than (4). Note also that 
both expressions do conserve probability, i.e. 



This is, of course, a very useful property because all charges carried by the 
valence quarks are then automatically conserved. 

Contrary to the leading twist expression (6), the distribution (7) exhibits 
transverse broadening as the density of the target increases (the scattered 
quarks are pushed out to larger q t ). Consider now the probability of in- 
elastic scattering (i.e. with color exchange) to small transverse momentum. 
This is given by expression (4,7), integrated from q t = to q t = A: 






(8) 



A 









a We repeat, however, that the transverse broadening of the distributions of released 
partons does play an important role since it destroys the coherence of the projectile 
wave function 7 and affects the fragmentation into hadrons. 
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Here, we have written only the leading term in A 2 /Q 2 , neglecting subleading 
power-corrections and exponentially suppressed contributions. Hence, soft 
forward inelastic scattering is power-suppressed in the black body limit 
because the typical transverse momentum is proportional to Qf. This 
steepens the longitudinal distribution dN/dxp of leading particles since 
partons with large relative momenta fragment independently 6 ' 7 . On the 
other hand, for low target density, the projectile's coherence is not destroyed 
completely and leading quarks may recombine, recovering the "leading- 
particle" effect observed in pp scattering at not too high energy. This 
recombination effect should be taken into account when modeling minimum 
bias pA collisions in order to ensure a smooth transition from the high- 
density to the low-density regime; our implementation is described and 
studied in more detail in section 3. 

Integrating over the transverse momentum of the scattered quark, the 
elastic and total scattering cross sections for quark- nucleus scattering are 13 : 



Clearly, when Q s /A — > oo, the cross section approaches the unitarity limit. 
2.2. Gluons 

Gluon bremsstrahlung dominates particle production at xf ^ 0.1. At very 
large transverse momentum, q t 3> Q s , the inclusive gluon distribution is 
given in collinear factorization by the usual gg — > gg LO hard scattering 
function convoluted with the DGLAP evolved leading-twist gluon distri- 
bution of the projectile and target. However, for the high-energy part of 
the air shower only the pt-integrated longitudinal distribution of hadrons 
matters (cf. section 4) which is dominated by fragmentation of gluons with 
transverse momenta up to <~ Qf. In that regime leading- twist perturbative 
QCD can not be applied reliably. 

Gluon radiation with transverse momentum q t ~ Qf in high-energy 
hadron- nucleus collisions has been discussed in detail in 19 ' 20 . The release of 
gluons from the hadronic wave functions can be described by convoluting 
the gluon distribution in the hadron with a (semi-)hard scattering cross 
section. The main qualitative features of the bremsstrahlung spectrum is 
that it flattens from ~ 1/qf for asymptotically large q t to l/g 2 for q t in 




(10) 



(11) 
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between the two saturation momenta; for q t — > 0, finally, it approaches a 
constant (up to logarithms) 21 . 

For the present purposes we require a simple ansatz that can be easily 
implemented in a Monte-Carlo model and, at the same time, does incorpo- 
rate the above features. A useful approach has been suggested in refs. 20 . 
The "fusion" of two gluon ladders gives rise to the bremsstrahlung spectrum 

<& 

E d%^ 4n Nf 1 ~[^ J dk t a s{k 2 t )Mxi,k 2 t )^ A {x 2 ,{q t ~k t ) 2 ) , (12) 

where 4>(x, Q 2 ) denotes the unintegrated gluon distribution function of the 
projectile hadron or target nucleus, respectively. It is related to the gluon 
density by 

Q 2 

xg(x,Q 2 ) = J dk 2 t <P(x,k 2 t ) . (13) 
Eq. (12) can be integrated by parts to read 

-^-2 = 47ra s ( g 2 ) jfrti ^S x i9h(xi,q 2 )x 2 g A (x 2 ,q 2 ) . (14) 
The ansatz from 20 for the infrared-finite gluon densities is 

x g(x, Q 2 ) cx — min(Q 2 , Q 2 (x)) (1 - x) 4 , (15) 

a s 

with a s evaluated at max(Q 2 ,Q 2 ). Note that for large Q 2 , the IE- 
dependence of the gluon distribution exhibits the conventional xg(x) <~ 
x~ A (l — x) 4 behavior; this follows from the evolution of the saturation mo- 
mentum 13 Q s (x) ~ l/x x with x. On the other hand, for small Q 2 and 
x, the above ansatz exhibits a slow logarithmic growth xg{x) <~ log x~ x 
only. In any case, we consider (15) to be a simple parameterization which 
exhibits some generic qualitative features expected from gluon saturation 
(e.g. that it is of order l/a s at small Q 2 and x) while, at the same time, it 
is roughly consistent with the DGLAP gluon distribution at large Q 2 and 
x. In fig. 1 we compare the parameterization (15), with the normalization 
constant fixed by the condition J dx xg(x, Q 2 ) = 0.5 and with Q 2 ~ x~ - 3 , 
to the CTEQ5 LO distribution 22 . 



b For fixed coupling evolution this is true for any x; for running coupling evolution it 
holds only for not too small x, see section 2.3. 
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Figure 1. Comparison of the gluon distribution from cq. (15) with CTEQ5 LO, in the 
DGLAP regime (for a proton). 

The constant of proportionality in eq. (15) was chosen in 20 such as to 
reproduce the overall normalization of the charged hadron rapidity distri- 
bution in d + Au collisions at BNL-RHIC energy. This was possible because 
the forward region was not considered. On the other hand, here we need 
to consider the entire solid angle (in momentum space) and, in particular, 
ensure conservation of the energy carried by the projectile. In our approach 
we therefore fix the overall number of radiated gluons by the condition of 
energy conservation. This is discussed in more detail in section 3 where we 
also show that the charged hadron multiplicity in the central region of pA 
collisions at BNL-RHIC energy agrees roughly with that from 20 and with 
available data. 

2.3. The saturation momentum as a function of impact 
parameter and rapidity 

The saturation momentum of the nucleus, Qf, must of course depend on 
the impact parameter as it basically measures the color charge density in 
the transverse plane. Hence, in the rest frame of the nucleus, the most naive 
estimate (neglecting shadowing 15 ) is that there are A times more valence 
quarks in a nucleus than in a nuclcon which are distributed over an area 
proportional to A 2 / 3 ; this then results in (Qf ) 2 ~ A 1 / 3 . More elaborate 
estimates lead to an additional factor equal to the logarithm of the mass 
number. 

For realistic nuclei with a non-uniform density distribution in the trans- 
verse plane we must replace, of course, the A 1 / 3 factor by the number of 
nuclcons from the target which interact with the projectile, Na- 



N A (b) = Aa in (s)T A (b) 



(16) 
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with Ta(1>) the nuclear profile function and <7i n (s) the energy dependent 
inelastic (non single-diffractive) cross section for the particular projectile 
species on protons. Fluctuations in A^i(6) are also taken into account in 
our Monte-Carlo implementation, as discussed in section 3. The impact 
parameter dependence of the nuclear saturation momentum is then taken 
to be 

Qf(b) = Ay/[1 + N A (b)} log(l + N A (b)) . (17) 

For a single nucleon, this corresponds to a saturation momentum on the 
order of A. It should be noted that at very high energies, deep in the black- 
body limit, the results depend only weakly on the above "initial condition" 
for Qf. However, when the saturation momentum is only moderately large 
(for example, for running coupling evolution, see below), the assumed de- 
pendence of Qf on b could play a role. Whether or not Q 2 S ~ Na is the 
most appropriate choice will be studied in more detail in the future. 

Next, we turn to the dependence of Q s on rapidity d y = \ogl/x. 
Eq. (17) provides the initial condition in (or near) the rest frame of the 
nucleus, y ~ 0, from valence quarks. As one moves away in rapidity phase 
space for gluon radiation opens up and so the gluon density grows; it is 
expected to saturate when it becomes of order l/a s 11 . For a recent review 
of evolution at small x see e.g. 23 . 

Model studies of deep inelastic scattering (DIS) on protons at HERA 
suggest 24 

Ql{x) ~ x~ x (18) 

with A w 0.3. This scaling relation can be obtained from the fixed coupling 
BFKL evolution equation for the scattering amplitude of a small dipolc. 
The BFKL equation is a linear QCD evolution equation which can not be 
applied in the high-density regime. Nevertheless, one can evolve the wave 
function of the target in rapidity y = logl/a; and ask when the dipole 
scattering amplitude becomes of order one, which leads to c 

Ql(y, b) = Q 2 s (yo, b) exp ca s y , (19) 

with a s = a s N c /-K and c w 4.84 a constant. Hence, LO fixed-coupling 
BFKL evolution predicts A' = ca s of order one, a few times larger than the 



c Normalized according to J d?bTA{b) = 1. 

d In this section, we measure the rapidity always relative to the parent hadron. 
° We write only the leading term proportional to y. 
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fit (18) to HERA phenomenology. A rcsummcd NLO BFKL analysis cor- 
rects this discrepancy and leads to A' much closer to the phenomenological 
value 25 . A similar observation is made in 5 where both log(l/a;) and logQ 
effects were considered. 

On the other hand one could also consider BFKL evolution with ad- 
hoc one-loop running of the coupling 234 : a s (Q 2 ) — 6 / \og(Ql(x) / Aq CD ) , 
which leads to 

Ql (V,b) = A 2 QCD exp ^2b c(y + y ) , (20) 

with 2& q/o = l°g 2 ( < 3o(^) 2 /AQ CD ). Insisting that (18) be valid at least in 
the y — > limit again provides us with a phenomenological value for the 
constant c in terms of the saturation momentum at y = 0. The form (20) 
leads to a notably slower growth of Q s at high energy. Specifically for 
central proton-nitrogen collisions at RHIC, LHC and GZK-cutoff energies 
(total rapidity y = 10.7, 17.3 and 26.0) the saturation momentum of the 
nucleus in the rest frame of the projectile hadron is Q s = 1.5, 5, 20 GeV for 
fixed coupling evolution, while for running coupling evolution it is Q s = 1, 
2.5,6 GeV, respectively. Clearly, cosmic ray interactions in our atmosphere 
should offer a realistic opportunity for distinguishing these scenarios. 



3. Monte-Carlo implementation 

We first generate a configuration of valence quarks according to the distri- 
bution (3,7), convoluted with the respective valence quark distribution of 
the projectile at the scale Qf 1 : 

Mx,Qf)C(q t ) (21) 



dxd 2 q t d 2 b 

with Qf is a function of both x and b. For this purpose, we employ the 
GRV94 parameterization of the parton distribution functions of a proton 26 
or a 7r +27 ; we assume isospin symmetry to deduce the distributions for other 
states. Also, as a rough approximation we take the valence quark distri- 
bution of the K + to be the same as that of the tt + , with the replacement 
d^ a. 

The remaining momentum is then used to generate a number of glu- 
ons according to the distribution (14). These gluons could be fragmented 
independently but this is not a good approximation when their transverse 



Again, we drop subleading terms that grow more slowly with rapidity. 
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momenta are soft. Rather, soft collinear gluons should be absorbed by the 
parent parton. 

The Lund string model 28 provides such an infrared safe fragmentation 
prescription. We order the produced gluons in rapidity and place them 
on strings between the valence quarks and the target nucleus, whose pre- 
cise configuration is not important. (Target fragmentation produces only 
low-energy particles which do not affect the properties of the air shower 
discussed here.) 

A baryon-nucleus collision produces three strings (two for a meson- 
nucleus collision) . However, when the invariant mass of any two of the three 
valence quarks is small, one cannot assume anymore that those strings frag- 
ment independently Rather, they will recombine to form a leading diquark, 
recovering the "leading particle effect" for low Qf. This effect should be 
taken into account in order to ensure a smooth transition from the regime 
of high target density (high energy, central collisions) to low target density 
(lower energy, large impact parameter). We model this by introducing a 
cut-off in invariant mass, 

m cut = to p = 0.77 GeV (22) 

below which two leading quarks are allowed to form a diquark. The effect 
on the ^-distribution of fast hadrons is shown below. 

Although soft gluon absorption and diquark recombination are taken 
into account in our Monte-Carlo implementation of scattering near the 
black body limit of QCD ("BBL"), it should nevertheless be clear that 
it is restricted to the high-density regime. For example, when Qf becomes 
small, the DGLAP leading-twist regime becomes important and one should 
use better approximations for the gluon densities than those from (15). 
Also, the fraction of diffractive and elastic events becomes sizable. 

A large amount of work has been done to develop models for this regime. 
S1BYLL 29 and QGSJET 30 , in particular, are commonly used to model air 
showers. We do not intend to duplicate those approaches here but rather to 
study whether anything could be learned about small- a; QCD from cosmic 
ray air showers. Hence, we couple our model to the standard pQCD leading 
twist event generator SIBYLL 2.1 such that the "BBL" Monte-Carlo treats 
the high density regime (large saturation momentum of the nucleus, i.e. high 
energy and/or small impact parameter) while SIBYLL handles peripheral 
or low energy collisions where the saturation momentum of the nucleus is 
not sufficiently large. 

It is clear, of course, that no sharp boundary between those regimes 
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exists and that this artificial separation is performed for purely technical 
reasons. It is therefore important to check that the results do not depend 
strongly on the precise location of the assumed boundary between low and 
high density. As already mentioned above, we do implement some effects 
into the BBL model which should facilitate a smooth transition to low 
densities, such as soft gluon absorption and diquark recombination. At the 
same time, the SIBYLL model also assists the transition to high target 
densities by implementing a low-p t cutoff for the DGLAP regime which 
grows rapidly with energy, see Engcl ct al 29 . We therefore expect that our 
results are not very sensitive to where exactly we perform the switch, as 
long as it occurs in a reasonable regime; this will be checked below. 

The saturation momentum of the nucleus provides an intrinsic scale 
for resolving the valence quark structure of the projectile, cf. eq. (21). In 
practice, this scale can not be too small because the Q 2 evolution of parton 
distributions in hadrons is normally obtained from DGLAP; standard PDF 
parameterizations typically require a minimal Q 2 on the order of 1 GeV 2 . In 
order not to distort the inclusive momentum distribution of valence quarks, 
we must therefore ensure that Qf (x) does not drop below this threshold at 
too large an x. Specifically, we require that the valence quark distribution 
be probed at least down to x = 10~ 3 : 

Q?(x = 1(T 3 , b) > Q min « 1 GeV . (23) 

In our Monte-Carlo approach, the collision is handled by either BBL or 
SIBYLL depending on whether this condition is met or not. The result- 
ing boundary between low and high density regimes appears reasonable; 
for example, for central collisions of protons on heavy targets like Au or 
Pb, the transition occurs just below BNL-RHIC energy, y/s = 200 GeV. 
On the other hand, minimum bias pp collisions essentially never pass the 
threshold (23), even at LHC energies and beyond. 

As a first check, we apply our model to RHIC energy which represents 
the highest presently available energy for proton- nucleus collisions. We 
compare the (pseudo-) rapidity distribution of inclusive charged hadrons to 
data by BRAHMS 31 . The fragmentation region of the target should be dis- 
regarded since no attempt has been made to treat that realistically. Given 
that no special tuning has been performed to fit these particular data s , we 
consider the qualitative agreement to be quite good. More importantly, we 



g For example regarding diquark recombination, string fragmentation or initial conditions 
for the evolution of the two saturation momenta 
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Figure 2. Left: Comparison of the BBL event generator (with running coupling evolu- 
tion) to RHIC data (E C M = 100 GcV) from d + Au collisions by the BRAHMS collab- 
oration 31 . We have scaled our results for p + Au (obtained with Na = 6 participants) 
by a factor of two. Right: our prediction for central p + Pb collisions at LHC energy, 
(assuming N A = 10 and E CM = 3000 GeV). 

note that both evolution scenarios (running and fixed coupling) can easily 
be made to fit the same data at this energy by somewhat readjusting the 
initial conditions for Q s (see for example 20 for a much better fit than ours 
with fixed-coupling evolution). Hence, RHIC energy is too low to reliably 
probe the evolution of Q s ; rather, results are mostly sensitive to the initial 
conditions. 

In the right panel we show our result for central p + Pb collisions at 
LHC energy, which roughly agrees with that from ref. 32 . Note, however, 
that ref. 32 considers the overall normalization to be a parameter, fixed from 
low-energy (RHIC) collisions, while in our approach it is determined auto- 
matically by momentum conservation. The similarity of dN/drj at central 
rapidities obtained via the two methods perhaps suggests that indeed the 
number of radiated gluons equals the maximum number allowed by kine- 
matics. 

Figure 3 shows the xf distribution of pions and nucleons for central 
proton- nitrogen events at 10 9 GeV. We plot on a logarithmic scale to show 
the effect at high xf- Very forward particle production is suppressed as 
compared to the pQCD model SIBYLL 2.1, which is a consequence of the 
break-up of the projectile into its partonic components. This behavior 
affects another key quantity for cosmic ray air showers, the so-called inelas- 
ticity, which is one minus the Feynman-x of the most energetic secondary 
hadron (shown below). Also note that in the high-density limit diquark 
recombination in the forward region is suppressed (see discussion below), 
and so the projectile proton mainly decays into a beam of leading mesons 7 . 
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Figure 3. Comparison of the SIBYLL 2.1 and BBL models for central p + 14 TV events at 
E = 10 9 GeV. For either model, we show the nucleon and pion spectra separately. One 
observes the suppression of forward nucleon production in the high-density limit, which 
is due to the complete breakup of the proton. 




It would be very useful to check this prediction in central p + Pb collisions 
at the LHC in order to confirm or rule out the basic mechanism of en- 
ergy degradation presented here, which is very important for cosmic ray air 
showers. 




10° 10' 
energy [GeV] 



10° 10° 
energy [GeV] 



Figure 4. Mean multiplicity of charged particles (left panel) and the inelasticity (right 
panel) as a function of lab energy, for minimum bias p + 14 N collisions. 



Fig. 4 compares the mean multiplicities of charged particles from the 
BBL model with fixed and running coupling evolution of the saturation 
scale, and the conventional models SIBYLL 2.1 and QGSJET01. Most signif- 
icant is the difference between fixed and running coupling evolution. For 
fixed coupling evolution Qf grows very large over a broad range of impact 
parameters (the radius of the black disc approaches the geometrical cross 
section of the target nucleus). Hence, even for minimum bias collisions 
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forward particle production is strongly suppressed in this case, and the re- 
maining energy is used for particle production at small xf- This explains 
the large multiplicities as compared to BBL with running coupling evolu- 
tion. On the other hand, for energies below 10 4 — 10 5 GeV there is little 
sensitivity to the evolution scenario for Qf and the results look rather sim- 
ilar. In SIBYLL 2.1, the growth of the multiplicity is "tamed" by a rapidly 
growing p t -cutoff for leading-twist hard processes. It should be kept in 
mind though that the multiplicity from large momentum transfer processes 
is power-law sensitive to the infrared cutoff. 




10 2 10 3 10 4 10 5 10 6 10 7 10 s 10 s 10 10 10 11 10 z 10 3 10 4 10 s 10 e 10 7 10 s 10 9 10'° 10 11 

Energy [GeV] Energy [GeV] 

Figure 5. The mean charged particle multiplicity and the inelasticity for minimum-bias 
p+ 14 N collisions using the combined BBL+SIBYLL 2.1 model for Q min = 0.7, 1.0 GeV. 

To check the sensitivity to the (artificial) boundary between low and 
high density from eq. (23), we compare results for the multiplicity and 
for the inelasticity for Q m i n = 0.7 GeV and Q m i n = 1 GeV in Fig. 5. A 
lower value for Q min leads to a higher fraction of BBL events, but these 
are then generated with lower Qf, and so are more similar to "soft" events 
from SIBYLL. In total, we see that there is little sensitivity of physical 
observablcs to the precise threshold between the models, as long as it is 
chosen within reasonable bounds. In the following we chose <5 m i n = 0.7 GeV 
as default for our calculations. 

In Fig. 6, finally, we show the effect of the diquark recombination mech- 
anism. We compare the production of protons and neutrons in central colli- 
sions to the case without recombination (m cut — 0). At relatively low ener- 
gies (E « 10 6 GeV) , one notices a suppression of forward baryon production 
when recombination is not taken into account, except for xf ~ 1- This very 
forward peak is in fact produced by elastic or diffractive events within the 
SIBYLL model, which still handles about 5% of all central collisions at this 
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Figure 6. Diquark recombination effect for p + 14 N collisions at various energies. 

energy if Q s obeys running coupling evolution 11 . At E = 10 8 GeV and 
above, essentially all central collisions occur near the black body limit, i.e. 
they pass criterion (23). The diquark recombination mechanism then al- 
lows more particles to be produced in the forward region since the momenta 
of the corresponding valence quarks are combined. Nevertheless, the effect 
is less important at higher energies since there Qf is already too high to 
allow for the production of a diquark system with low invariant mass. 

4. Air showers 

In section 4.1 we give a brief introduction into concepts and observables in 
cosmic ray physics for readers from other fields. More detailed discussions 
can be found in dedicated textbooks such as the book by Gaisser 33 . 

In section 4.2 we discuss general aspects of air shower simulation and, in 
particular, present the so-called cascade equations 34 employed here to solve 
for the longitudinal shower profile. From those equations, we can cleanly 
identify which "input" is required from QCD for their solution and, indeed, 
which properties of high-energy hadronic interactions actually influence the 
characteristics of very high energy air showers. 

4.1. Introduction to air showers 

Due to a very low flux at high energies, cosmic rays are detected indirectly 
by the measurement of air showers. These are cascades of particles produced 



h We note that the actual configuration of nuclcons in the target is generated randomly 
in each event by SIBYLL according to the appropriate nuclear density profile. Hence, 
fluctuations in the number of target participants are taken into account. 
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Figure 7. Longitudinal profiles and lateral distribution function of typical showers; here 
E = 10 10 GeV. 



by the interaction of the primary cosmic ray and subsequent secondaries 
with air nuclei. An air shower can be structured into 3 parts: a hadronic, 
an electromagnetic and a muonic part. Hadrons are produced in collisions 
with air-nuclei. Most of the electromagnetic part is induced by 7r°-decays, 
which have a short life-time and decay instantly up to 10 19 eV. Muons are 
produced by decays of charged pions and kaons, but since their decay-length 
is much longer only low energies particles decay while at higher energies 
collision with air nuclei dominates. Once produced, muons propagate with 
little interaction (mostly energy-loss) through the atmosphere until they 
decay or reach the ground. The most prominent fraction of a shower is 
the electromagnetic part. A rule of thumb is that the number of electrons 
and positrons at the maximum of the shower is approximately 60% of the 
primary energy Eq measured in GeV (a 10 11 GeV shower produces about 
60 billion particles !). 

There are two basic observables associated with air showers, which are 
measured by experiments: the longitudinal shower profile and the lateral 
distribution functions. The longitudinal profile is the number of charged 
particles measured along the shower-axis. One typically expresses this as a 
function of slant-depth, which is the density of the atmosphere integrated 
along the shower-axis: X = PAir(l)dl. For a vertical shower, X ranges 
from zero (top of the atmosphere) to 1020 g/cm 2 (sea level). For inclined 
showers with polar angle up to 60°, the slant depth is related to the vertical 
depth by X = X v / cos(0); at larger angles one has to take into account the 
curvature of the earth. Typical shower profiles are shown in Fig. 7. The 
position where the profile reaches the maximum is defined as the shower 
maximum A max , the number of charged particles is called the shower size 
N max . For a fixed energy these values fluctuate quite substantially, which is 
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due to the fact that the depth of the first collision can vary according to the 
cross section. One therefore usually compares the mean A max for a given 
primary and energy Eq. An important observation is that A max oc log(-Eo) 
and A max oc Eq- To first approximation, the shower of a nucleus can be 
considered to be a superposition of A independent nucleon-initiated show- 
ers, each carrying an energy Eq /A {Eq is the total energy of the primary, 
not per nucleon). The mean A^ ax of an iron- induced shower is therefore 
lower than that of a proton induced shower, roughly corresponding to the 
mean A max of a proton shower at energy Eq /A. Experiments measure the 
longitudinal profile via the emitted fluorescence light of nitrogen as the 
shower swipes through the atmosphere. The energy of the primary cosmic 
ray is then proportional to the total number of charged particles, which is 
determined by integrating the profile. 

A very important quantity for A max is the inelastic cross section of 
a particle on air. It determines the mean free path in the atmosphere. A 
significant amount of uncertainty in models for the longitudinal distribution 
of air showers is due to this variable. 

The other observable, the lateral distribution function (LDF), describes 
the density of particles measured on the ground as a function of the distance 
from the shower axis (hence in the shower plane) for given particle types 1 . 
Most of the lateral spread is generated by low-energy scattering of the 
electromagnetic part of an air shower. Hadrons do not spread out very 
much to large distances, only the low energy ones influence the tail of the 
LDF by producing tt° at large angles or distances j . Typical LDFs are 
depicted in Fig. 7 (right panel). They follow approximately a power law. 
Just as A max , the slope of the LDF also fluctuates. Showers induced higher 
in the atmosphere lead to flatter LDFs since they spread out over a larger 
radial distance. Empirically, one finds that these fluctuations cancel at some 
distance from the shower axis. Experiments exploit this property to extract 
the primary energy of the cosmic ray, which is taken to be proportional to 
the density at some distance from the axis. The proportionality constant 
is normally computed from simulations. 

When studying high energy particle physics with air showers it is im- 
portant to notice that the high-energy part of the shower (i.e. the first 



'Experiments measure the density in terms of response of whatever detector they use, 
e.g. scintillation or Cerenkov light, and normalize by the average signal of atmospheric 
muons, which are used for calibration. 

jflcncc, LDFs constrain mostly the low-energy hadronic interaction models 35 . 
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collisions) is almost purely longitudinal, which follows from simple kine- 
matics. The targets (air-nuclei) are at rest and the projectiles have huge 
7-factors, resulting in very small scattering angles. Furthermore, forward 
scattering is most important in air showers, since large-x^ particles carry 
most of the energy. This implies, for example, that high-pi QCD jets at 
mid-rapidity do not influence the longitudinal shower profile substantially 
nor do they contribute significantly to the lateral spread (see also ref. 36 ). 

4.2. Simulation of air showers 

The simulation of air showers is crucial for cosmic ray physics, since they 
are needed for the interpretation of experimental data. Given a hadronic 
interaction model k , and models for electromagnetic and muonic interactions 
one could just follow each particle and subsequent secondaries individually. 
Of course, since N cx Eq, this would require huge amounts of computing 
time at very high energies. Therefore, Hillas introduced the thinning algo- 
rithm: below a given energy threshold, i.e. for E < f t h x Eq, only one single 
secondary particle from a collision is followed, but it is attributed a higher 
weight. However, for a large thinning level fth, this introduces artificial 
fluctuations into the air shower. 

On the other hand, the fact that at high energies the lateral part 
of the hadronic shower can be neglected suggests another efficient ap- 
proach to solve this problem, which is based on one-dimensional transport 
equations 34 : 



Here, h n (E, X)dE is the number of particles of type n at altitude X in the 
given energy range [E, E + dE]; the functions W mn (E', E) are the energy- 
spectra dN/dE of secondary particles of type n in a collision of hadron m 
with air; D mn (E', E) are the corresponding decay functions; d n = m n / (ct„) 
is the decay constant, and X n (E) cx l/Vinei is the mean free path of the par- 
ticle. The first term in (24) with the minus sign accounts for particles dis- 
appearing by either collisions or decays, whereas the source term accounts 



dh n {E,X) 
dX 



-h n (E,X) 



X n (E) Ep(X) 



(24) 




m 



k Typically a Monte-Carlo event generator which generates complete final states and 
accounts for fluctuations. 
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for production of secondary particles by collisions or decays of particles at 
higher energies. Primary particles appearing in eq. (24) are nucleons (pro- 
tons, neutrons and their anti-particles), charged pions, charged and neutral 
kaons. In addition, we have as secondaries n°s and photons (as direct de- 
cay products from rj mesons, for example) which feed the electromagnetic 
cascade and muons as decay products of charged mesons. 

To summarize, the basic ingredients for constructing longitudinal pro- 
files of air showers are the inclusive spectra dN n /d,XF of the non-strongly 
decaying particles and their inelastic cross sections, which determine the 
mean free path. The electromagnetic cascade can be treated in a similar 
way 5 ' . 

The first few interactions in an air shower are the main source of fluctu- 
ations in X maK (and, accordingly, in the LDF). Since the cascade equations 
cannot account for those (they solve for a mean shower) one could treat 
the high energy part by a traditional Monte-Carlo method. This is the so- 
called hybrid approach to air shower simulations. On the other hand, if one 
solves the cascade equations without fluctuations 38 one can still reproduce 
the average X max to within a few g/cm 2 . 



4.3. Sensitivity of X max to the xp distribution 

Finally, we analyze which region of the ^-distribution is most important 
for the mean X max of an air shower. From the simple argument that forward 
particles carry most of the energy it should be clear that the high x p region 
is important. Our goal here is to quantify this statement somewhat. 

Given dN n /dxp distributions of secondaries, to study the sensitivity 
of X max to various regions of xp we solve the cascade equations with a 
modified distribution: 

dN n dN n 

► - — (1 + e) for x F < x F . (25) 

dx F dxp 

That is, we enhance or suppress the spectra at sj? < x° F relative to the 
default reference distributions, depending on the sign of e. At the same 
time, we suppress or enhance particles at xp > x° F in such a way as to 
conserve the total energy: 



(26) 
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with 

x° 

E fdx F E n (x F ) fa 
e ' = e JL^ . (27 ) 

E / 

n -r-0 

Note that the dN n /dxF are, of course, energy dependent while we take the 
1 + e factor in (25) to be constant. Also, we do not modify the inelastic cross 
section (i.e. the mean free path in the atmosphere), just the x^-distribution 
of secondaries in an inelastic event. We then solve eqs. (24) to determine the 
change of X max relative to that for the reference distributions as a function 
of both e and x° F . 




Figure 8. The shift of -X" max as a function of e for various x p . This shows that the 
shower maximum is sensitive to mainly the forward region, xp > 10~ 3 . 



The result is shown in Fig. 8, assuming a proton primary with energy 
Eq = 10 10 GcV. The reference dN n /dxp distributions were taken from 
QGSJET01. We observe that AX maK is approximately linear in e. For large 
x F , for example =0.1, there is a significant shift AX max w 140e. For e < 
we suppress the Xf < 0.1 region and, by energy conservation, enhance the 
large-XF part; this leads to deeper penetration into the atmosphere, i.e. to 
larger X max . In turn, suppression of forward particle production (e > and 
e' < 0) leads to decreasing X maK . 

However, one can also observe that X max becomes independent of e for 
x F < 10 -3 . This shows that the small-x^ part of the distribution has no 
influence on the shower maximum (for fixed cross section). For comparison, 
we note that a particle produced at mid-rapidity (in a collision with energy 
10 10 GeV) with an energy of about the proton mass has x F ~ 10~ 5 . 
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5. Application of the BBL to air showers 
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Figure 9. Left panel: mean X max (_Eo) for p and Fe induced showers from SIBYLL 2.1, 
and for p primaries with BBL for fixed and running coupling evolution; data are from 
HIRES 2 ; figure from 8 . Right panel: mean lateral distribution function scaled by the 
LDF found by AGASA. 



To apply our interaction model to air showers, we tabulated dN n /dE 
distributions of the primary particles appearing in an air shower. We then 
employ the SENECA 34 model to solve the cascade equations (24). The 
hadron-air inelastic cross sections are taken as parameterized in SIBYLL 2.1. 
In Fig. 9 we compare the results for fixed and running coupling evolution to 
those obtained with the SIBYLL 2.1 model. One notices a huge difference 
between fixed and running coupling evolution scenarios. The saturation 
momentum in the former case is so high that forward scattering is very 
strongly suppressed over a broad range of impact parameters; also, the 
discrepancy between those evolution scenarios at the highest energies is 
strongly amplified by subsequent hadronic collisions in the cosmic ray cas- 
cade. Consequently, for fixed coupling evolution the shower is absorbed very 
early in the atmosphere. Hence, if we assume a hadronic primary, then the 
HIRES 2 data excludes this scenario, since it would require hadrons lighter 
than protons. This is a novel result as present accelerator data could not 
rule out such a rapid growth of the gluon density (see e.g. 24,20 ); it illustrates 
the ability of cosmic ray air showers to provide observational constraints 
on small- a; QCD 8 . 

The running coupling result, on the other hand, is compatible with 
those data and with a light composition. The results from this model are 
similar to those from SIBYLL 2.1 or QGSJET01 (for proton primaries) within 
present theoretical uncertainties. Nevertheless, our results show that the 
effects discussed here make near-GZK proton-induced air showers look more 
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similar to nucleus-induced cascades from leading-twist models and so favor 
a lighter composition near the cutoff. 

Lateral distribution functions obtained with both evolution scenarios 
are shown in the right panel of Fig. 9. The LDFs were computed for the 
AGASA experiment and include a full detector response simulation of the 
plastic scintillators. The results are scaled by the empirical formula which 
describes the data up to highest energies quite well: 

™-^)>£)">te)T- ,28) 

with a — 1.2, 5 = 0.6, Rm — 91.6 m and r\ — 3.84 for vertical showers 39 . 
The parameter C is adjusted using the energy conversion formula (13) from 
Rcf. 39 , 

E= 2.17 x 10 17 5(600 m)eV, (29) 

which is valid for the average altitude of the AGASA array, =667 m. The 
comparison in Fig. 9 shows that the LDF obtained for fixed coupling evolu- 
tion is much flatter than that for running coupling evolution, which in turn 
agrees better with the data (notice that in the figure, the theoretical curves 
are scaled by the data). This is consistent with our finding for A max , as 
discussed above. When the shower is absorbed earlier in the atmosphere 
then it spreads out to larger radial distances from the shower axis. 

6. Conclusion and Outlook 

In this paper we pointed out that atmospheric air showers induced by the 
highest energy cosmic rays are sensitive to QCD interactions at extremely 
small momentum fractions x where nonlinear effects arc expected to play a 
major role and lead to unitarization of partonic scattering cross sections. In 
turn, this means that cosmic rays air showers can provide valuable insight 
and observational constraints for the strong-field regime of QCD. As an 
example, we have shown that present data on A max (i?) already exclude that 
the rapid ~ 1/x 3 growth of the saturation boundary (which is compatible 
with RHIC and HERA data) persists up to GZK cutoff energies 8 . 

The model used here for quantitative calculations can be improved in 
many ways, for example by incorporating more advanced estimates for the 
small- x gluon densities obtained from the approaches of refs 4 ' 5 . From the 
point of view of learning about cosmic rays from small- a; QCD it could 
be interesting to extend the studies to nucleus- nucleus collisions and to 
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perform a composition analysis near the cutoff. This might also be relevant 
for physics well below the cutoff, in the region above the "knee" 40 (E ~ 
5 • 10 6 GeV) because for nuclei nonlinearities should set in at lower energies 
already. 

Cosmic ray air showers offer several important advantages over labo- 
ratory experiments: first of all, of course, their energies can exceed those 
of accelerators (even LHC) by far. Second, many properties of extensive 
air showers are sensitive mainly to the forward region and to transverse 
momenta about (pt), which means that they probe extremely small x in 
the target nucleus. Finally, an air shower develops via several subsequent 
collisions and so any "distortion" of the momentum-space distribution of 
secondaries from high-density effects is strongly amplified (essentially raised 
to the power of the number of collisions) . 

On the other hand, unlike air shower detectors accelerator experiments 
can control key parameters of the interaction. For example, aside from 
collision energy and centrality one can also chose various projectiles and 
targets, from protons over light nuclei up to very massive nuclei such as 
gold or lead. Central collisions on lead at LHC energy should provide 
similar gluon densities as those on air at cutoff energies. Hence, fruitful 
lessons regarding small- a; QCD will hopefully emerge from both cosmic ray 
and accelerator data in the future. We emphasize that crucial data to be 
obtained at the LHC is not limited to the total proton-proton cross section 
but includes xf distributions of secondaries fromp+A collisions in both the 
central and forward regions. The latter would allow us to study the energy 
degradation mechanism in central collisions, which plays an important role 
for cosmic ray air showers. 
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